Symcap
Teegan Innis
September 21, 2016
Project Summary
This project aimed to look at the spatial variability of Symbiodinium clades C and D in the Kane’ohe Bay, O’ahu, Hawai’i population of Montipora capitata. We investigated the distribution of symbionts at scales ranging from location in the bay to location on an individual reef. We also looked at differences among reef types (fringing vs. patch), colony color morph (brown vs. orange) and depth. Heterogeneous mixtures of symbiont clades were considered in the analysis for spatial patterns. By investigating spatial variability of Symbiodinium, we furthered the understanding of stress-response potential in Kane’ohe Bay.
Library Packages
setwd("~/Symcap")
library(data.table)
library(lsmeans)
library(devtools)
library(plyr)
library(reshape2)
library(popbio)
library(RgoogleMaps)
library(plotrix)
library(zoo)
library(rgdal)
library(car)
library(scales)
library(png)
library(pixmap)
library(ecodist)
library(cluster)
library(fpc)
library(clustsig)
library(foreign)
library(nnet)
library(ggplot2)
library(mlogit)Import/Merge Field and qPCR Data
Coral_Data <- read.csv("Coral_Collection.csv")
Coral_Data$Depth..m. <- as.numeric(as.character(Coral_Data$Depth..m.))
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path = "qPCR_data", pattern = "txt$", full.names = T)
Mcap <- steponeR(files = Mcap.plates, delim="\t",
target.ratios=c("C.D"),
fluor.norm = list(C=2.26827, D=0),
copy.number=list(C=33, D=3))
Mcap <- Mcap$result
Mcap <- Mcap[grep("+", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("NTC", Mcap$Sample.Name, fixed = T, invert = T), ]
Mcap <- Mcap[grep("PCT", Mcap$Sample.Name, fixed = T, invert = T), ]
colnames(Mcap)[which(colnames(Mcap)=="Sample.Name")] <- "Colony"
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]
Mcap$C.D[which(Mcap$C.reps<2)] <- -Inf
Mcap$C.D[which(Mcap$D.reps<2)] <- Inf
Mcap <- Mcap[with(Mcap, order(Colony)), ]
Mcap$propC <- Mcap$C.D / (Mcap$C.D+1)
Mcap$propD <- 1-Mcap$propC
Mcap$propD[which(Mcap$C.D==-Inf)] <-1
Mcap$propC[which(Mcap$C.D==-Inf)] <-0
Mcap$propD[which(Mcap$C.D==Inf)] <-0
Mcap$propC[which(Mcap$C.D==Inf)] <-1
Mcap$Dom <- ifelse(Mcap$propC>Mcap$propD, "C", "D")
Symcap<-merge(Coral_Data, Mcap, by="Colony", all=T)
Symcap <- Symcap[with(Symcap, order(Colony)), ]
Symcap$Mix <- factor(ifelse(Symcap$propC>Symcap$propD, ifelse(Symcap$propD!=0, "CD", "C"), ifelse(Symcap$propD>Symcap$propC, ifelse(Symcap$propC!=0, "DC", "D"), NA)), levels = c("C", "CD", "DC", "D"))
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")Adjust Depth by Mean Sea Level
To account for the influence of tides, depth was adjusted according to the difference in sea level from the mean sea level on each collection day at 6-minute intervals. Mean sea level was obtained from NOAA tide tables for Moku o Lo’e.
JuneTide=read.csv("Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("Station_1612480_tide_ht_20160701-20160731.csv")
AugustTide=read.csv("Station_1612480_tide_ht_20160801-20160812.csv")
Tide<-rbind(JuneTide, JulyTide, AugustTide)
Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
Symcap$Time2 <- as.POSIXct(paste(as.character(Symcap$Date), as.character(Symcap$Time)),
format="%m/%d/%y %H:%M", tz="Pacific/Honolulu")
Symcap$Time=Symcap$Time2
# Add estimated times for missing values
Symcap$Time[170] <- as.POSIXct("2016-06-14 12:07:00")
Symcap$Time[177] <- as.POSIXct("2016-06-14 12:20:00")
Symcap$Time[178] <- as.POSIXct("2016-06-14 12:22:00")
Symcap$Time[180] <- as.POSIXct("2016-06-14 13:08:00")
Symcap$Time[187] <- as.POSIXct("2016-06-14 12:42:00")
Symcap$Time[188] <- as.POSIXct("2016-06-14 12:44:00")
Symcap$Time[206] <- as.POSIXct("2016-06-16 13:10:00")
Symcap$Time[208] <- as.POSIXct("2016-06-16 13:24:00")
Symcap$Time[211] <- as.POSIXct("2016-06-16 12:37:00")
Symcap$Time[218] <- as.POSIXct("2016-06-16 12:27:00")
Symcap$Time[448] <- as.POSIXct("2016-07-16 13:32:00")
Round6 <- function (times) {
x <- as.POSIXlt(times)
mins <- x$min
mult <- mins %/% 6
remain <- mins %% 6
if(remain > 3L) {
mult <- mult + 1
} else {
x$min <- 6 * mult
}
x <- as.POSIXct(x)
x
}
Symcap$Time.r <- Round6(Symcap$Time)
Tide$Time.r <- Tide$Time
merged<-merge(Symcap, Tide, by="Time.r", all.x=T)
merged$newDepth <- merged$Depth..m.- merged$TideHTCoral Collection Map
KB <- c(21.46087401, -157.809907)
KBMap <- GetMap(center = KB, zoom = 13, maptype = "satellite", SCALE = 2, GRAYSCALE = FALSE)
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude, col=153, pch=21, bg="#FF7F50", lwd=2)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
par(new=T, mar=c(11.8,19,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()
Chi-Squared
Significant Differences
Dominant Symbiont by Reef Area
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Dom, Symcap$Reef.Area)
chisq.test(results)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 136.26, df = 1, p-value < 2.2e-16prop.table(results, margin = 2)##
## Slope Top
## C 0.7767857 0.3294574
## D 0.2232143 0.6705426par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Symbmiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)
Dominant Symbiont by Color Morph
results=table(Symcap$Dom, Symcap$Color.Morph)
chisq.test(results)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 164.96, df = 1, p-value < 2.2e-16prop.table(results, margin = 2)##
## Brown Orange
## C 0.8896321 0.4103194
## D 0.1103679 0.5896806par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Color Morph", ylab = "Symbiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)
Symbiont Community Composition by Color Morph
results=table(Symcap$Mix, Symcap$Color.Morph)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 167.44, df = 3, p-value < 2.2e-16prop.table(results, margin = 2)##
## Brown Orange
## C 0.762541806 0.361179361
## CD 0.127090301 0.049140049
## DC 0.107023411 0.570024570
## D 0.003344482 0.019656020par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Color Morph", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)
Symbiont Community Composition by Reef Area
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Mix, Symcap$Reef.Area)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 138.97, df = 3, p-value < 2.2e-16prop.table(results, margin = 2)##
## Slope Top
## C 0.678571429 0.275193798
## CD 0.098214286 0.054263566
## DC 0.214285714 0.651162791
## D 0.008928571 0.019379845par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Reef Area", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)
Symbiont Community Composition by Dominant Symbiont
results=table(Symcap$Mix, Symcap$Dom)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 707, df = 3, p-value < 2.2e-16prop.table(results, margin = 2)##
## C D
## C 0.86635945 0.00000000
## CD 0.13364055 0.00000000
## DC 0.00000000 0.96703297
## D 0.00000000 0.03296703par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray 10", "gray 85", "gray 40", "gray100"), xlab = "Dominant Symbiont", ylab = "Symbiont Mixture Proportion")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray85", "gray40", "gray100"), inset = c(-.2, 0), xpd = NA)
When looking at D-only colonies, the Ct values are on the higher end of the spectrum (35 or greater) indicating poor amplification and the potential for C amplification being pushed back in the cycle. This is supported by the fact that 5 of the 9 D-only colonies had C present in 1 qPCR replicate.
Proportion of D When Present in Mixture
propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]
hist(propD, xlab = "Proportion of Clade D", ylab = "Number of Samples", main = "", col = "gray75")
propDHist <- subset(merged, propD > 0 & propD < 1)
propDHist$propD <- cut(propDHist$propD, breaks = 5)
DCol=table(propDHist$Color.Morph, propDHist$propD)
z <- with(subset(merged, propD==0), table(Color.Morph))
o <- with(subset(merged, propD==1), table(Color.Morph))
DCol <- cbind(z, DCol, o)
par(mar=c(2, 4, 2, 0))
bars <- barplot(DCol, xlab = "Clade D Proportion", ylab = "Number of Samples",
main = "", col = c(alpha("sienna", 0.55), alpha("orange", 0.55)),
axisnames = FALSE, space = c(0,0.3,0,0,0,0,0.3))
axis(side = 1, at = c(1.3, 2.3, 3.3, 4.3, 5.3, 6.3), labels = c(">0.0", 0.2, 0.4, 0.6, 0.8, "<1.0"))
axis(side = 1, at = c(0.5, 7.2), labels = c("All C", "All D"), tick = FALSE)
Color Morph by Reef Area
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Color.Morph, Symcap$Reef.Area)
chisq.test(results)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 81.109, df = 1, p-value < 2.2e-16prop.table(results, margin = 2)##
## Slope Top
## Brown 0.5523385 0.2015504
## Orange 0.4476615 0.7984496par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Color Morph Proportion")
legend("topright", legend=c("Brown", "Orange"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)
Color Morph by Reef ID
results=table(Symcap$Color.Morph, Symcap$Reef.ID)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 50.691, df = 24, p-value = 0.001156prop.table(results, margin = 2)##
## 10 13 14 18 19 20
## Brown 0.4000000 0.4000000 0.5333333 0.2333333 0.3666667 0.6285714
## Orange 0.6000000 0.6000000 0.4666667 0.7666667 0.6333333 0.3714286
##
## 21 25 26 30 38 42
## Brown 0.4333333 0.6000000 0.4242424 0.1333333 0.2857143 0.5142857
## Orange 0.5666667 0.4000000 0.5757576 0.8666667 0.7142857 0.4857143
##
## 46 5 Deep F1-46 F2-25 F3-18
## Brown 0.5000000 0.3666667 0.5172414 0.4000000 0.6000000 0.2500000
## Orange 0.5000000 0.6333333 0.4827586 0.6000000 0.4000000 0.7500000
##
## F4-34 F5-34 F6-Lilipuna F7-18 F8-10 F9-5
## Brown 0.4400000 0.2500000 0.5000000 0.2000000 0.4000000 0.4000000
## Orange 0.5600000 0.7500000 0.5000000 0.8000000 0.6000000 0.6000000
##
## HIMB
## Brown 0.6000000
## Orange 0.4000000Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY) <- XY$Reef.ID
XY <- XY[, -1]
XY <- na.omit(XY)
apply(XY, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["Orange"], reef["Brown"]), radius = 7, col = c("orange", "sienna"))
})
Mantel Test for Spatial Autocorrelation
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
col.dists <- dist(XY$Brown)
set.seed(12456)
mantel(col.dists~reef.dists)## mantelr pval1 pval2 pval3 llim.2.5% ulim.97.5%
## -0.09132955 0.87800000 0.12300000 0.24500000 -0.14828356 -0.02809549Non-Significant Differences
Dominant Symbiont by Reef Type
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 0.77593, df = 1, p-value = 0.3784Dominant Symbiont by Bay Area
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 1.8678, df = 2, p-value = 0.393Dominant Symbiont by Reef ID
results=table(Symcap$Dom, Symcap$Reef.ID)
chisq.test(results)KB <- c(21.46087401, -157.809907)
KBMap <- GetMap(center = KB, zoom = 13, maptype = "satellite", SCALE = 2, GRAYSCALE = FALSE)
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propDom=table(Symcap$Dom, Symcap$Reef.ID)
propDom=prop.table(propDom, margin = 2)
propDom <- as.data.frame.matrix(propDom)
props <- data.frame(t(propDom))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY) <- XY$Reef.ID
XY <- XY[, -1]
XY <- na.omit(XY)
apply(XY, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["C"], reef["D"]), radius = 7, col = c("blue", "red"))
})
Mantel Test for Spatial Autocorrelation
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propDom=table(Symcap$Dom, Symcap$Reef.ID)
propDom=prop.table(propDom, margin = 2)
propDom <- as.data.frame.matrix(propDom)
props <- data.frame(t(propDom))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
dom.dists <- dist(XY$C)
set.seed(12456)
mantel(dom.dists~reef.dists)## mantelr pval1 pval2 pval3 llim.2.5% ulim.97.5%
## 0.15974663 0.04000000 0.96100000 0.04500000 0.08875833 0.24328773pval1 indicates that there is a positive correlation between distance and difference among reefs in terms of dominant symbiont clade. Reefs that are closer in proximity to each other are more similar in terms of symbiont dominance. A positive mantelr value indicates a positive correlation.
Clustering of Similar Reefs
doms <- hclust(dom.dists)
plot(doms, labels = XY$Reef.ID)
a <- hclust(reef.dists*dom.dists)
plot(a, labels = XY$Reef.ID)
domsk <- pamk(data = dom.dists, critout = TRUE)
domsk
km <- kmeans(dom.dists, centers = 5)
km
df <- cbind(XY, km$cluster)
km$cluster
PlotOnStaticMap(KBMap, df$Latitude, df$Longitude, col="red", pch=as.character(km$cluster), lwd = 2)Dominant Symbiont Plots
Raw Values
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["raw.c"], reef["raw.d"]), radius = 7,
col = c("red", "blue"))
})
Depth Adjustment
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["c.adj"], reef["d.adj"]), radius = 7,
col = c("red", "blue"))
})
Depth Interaction Adjustment
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["c.adj.int"], reef["d.adj.int"]), radius = 7,
col = c("red", "blue"))
})
Depth Adjustment for Color Morph and Dominant Symbiont
A multinomial logistic regression was performed on the interaction of color morph and dominant symbiont to discount the influence of depth on spatial distribution.
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID),
FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c
manteltable = table(merged$Dom, merged$Color.Morph, merged$Reef.ID)
nc <- aggregate(interaction(merged$Color.Morph, merged$Dom),
by=list(merged$Reef.ID), FUN=table)
nc <- data.frame(Reef.ID=as.character(nc[,1]), prop.table(nc[,2], margin=1))
XY3 <- merge(XY2, nc, by="Reef.ID", all=T)
merged$ColDom <- interaction(merged$Color.Morph, merged$Dom)
model <- multinom(ColDom~newDepth+Reef.ID, merged)
means <- lsmeans(model, specs = "Reef.ID")
means
pp <- fitted(model)
newdat <- data.frame(Reef.ID = levels(merged$Reef.ID),
newDepth = mean(merged$newDepth, na.rm=T))
pred <- predict(model, newdata = newdat, "probs")
new <- data.frame(Reef.ID=as.character(newdat[,1]), pred)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["Brown.C.x"], reef["Orange.C.x"], reef["Brown.D.x"],
reef["Orange.D.x"]), radius = 7,
col = c("turquoise4", "steelblue1", "tomato", "coral"))
})
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("turquoise4", "steelblue1", "tomato", "coral"))
par(new=T, mar=c(15,22,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()
Raw proportion values for Color Morph and Dominant Symbiont at each reef not discounting the effect of depth.
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["Brown.C.y"], reef["Orange.C.y"], reef["Brown.D.y"],
reef["Orange.D.y"]), radius = 7,
col = c("turquoise4", "steelblue1", "tomato", "coral"))
})
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("turquoise4", "steelblue1", "tomato", "coral"))
par(new=T, mar=c(15,22,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()
Adjusted values of Dominant Symbiont and Color Morph adjusted for the influence of depth.
Mantel Test for Spatial Autocorrelation
reef.dists <- dist(cbind(XY4$Longitude, XY4$Latitude))
dom.dists <- bcdist(cbind(XY4$Brown.C.y, XY4$Orange.C.y, XY4$Brown.D.y, XY4$Orange.D.y))
set.seed(12456)
mantel(dom.dists~reef.dists)## mantelr pval1 pval2 pval3 llim.2.5%
## 0.049858793 0.227000000 0.774000000 0.483000000 -0.008622074
## ulim.97.5%
## 0.111996212Color and Dominant Symbiont by Bay Area
When considering the effect of bay area on the interaction of color morph and dominant symbiont, the proportion of Brown colonies dominated by D increases as you move from the north to the south of the bay. The proportion increases almost 4x, yet the small number of colonies makes this non-compelling.
Type=table(merged$ColDom, merged$Bay.Area)
prop.table(Type)##
## Central Northern Southern
## Brown.C 0.124645892 0.113314448 0.138810198
## Orange.C 0.063739377 0.096317280 0.076487252
## Brown.D 0.014164306 0.007082153 0.025495751
## Orange.D 0.116147309 0.106232295 0.117563739chisq.test(Type)##
## Pearson's Chi-squared test
##
## data: Type
## X-squared = 12.95, df = 6, p-value = 0.04384Color Morph by Reef Type
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 1.1302, df = 1, p-value = 0.2877Color Morph by Bay Area
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 3.964, df = 2, p-value = 0.1378Symbiont Community Composition by Reef Type
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 2.8371, df = 3, p-value = 0.4174Symbiont Community Composition by Bay Area
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 8.9692, df = 6, p-value = 0.1753Symbiont Community Composition by Reef ID
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 87.127, df = 72, p-value = 0.1081Logistic Regression
Significant Effect
Dominant Symbiont by Depth
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
Dom1 <- subset(merged, !is.na(newDepth) & !is.na(Dominant))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 705 942.15
## newDepth 1 129.04 704 813.10 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1logi.hist.plot(Dom1$newDepth, Dom1$Dominant, boxp = FALSE, type = "hist", col="gray", xlabel = "Depth (m)", ylabel = "", ylabel2 = "")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "C D", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of clade C Symbiont", line = 3, cex = 1)
merged$DepthInt <- cut(merged$newDepth, breaks = 0:13)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results##
## (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
## 0 148 33 28 22 14 4 3 2 2 1 0
## 1 73 52 84 91 47 26 15 12 16 8 5
##
## (11,12] (12,13]
## 0 1 0
## 1 0 1props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.2, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Dominant Symbiont Proportion")
Bars indicate relative frequency of clade C vs. D dominance at 1m depth intervals when pooling all samples collected.
Mixture by Depth
merged$Mixture <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 1, 0)
merged$Mixture2 <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 0, 1)
results=glm(Mixture~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Mixture
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 705 973.27
## newDepth 1 93.481 704 879.79 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results=table(merged$Mixture2, merged$DepthInt)
results##
## (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
## 0 153 44 40 27 18 8 4 2 5 2 1
## 1 68 41 72 86 43 22 14 12 13 7 4
##
## (11,12] (12,13]
## 0 1 0
## 1 0 1props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Mix", "Non-Mix"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.23, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Mixture~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Mixture Proportion")
Bars indicate relative frequency of Mixture vs. Non-Mixtures at 1m depth intervals when pooling all samples collected.
Color Morph by Depth
merged$Color <- ifelse(merged$Color.Morph=="Orange", 1, 0)
results=glm(Color~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Color
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 706 963.85
## newDepth 1 78.717 705 885.14 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Color <- subset(merged, !is.na(newDepth) & !is.na(Color))
logi.hist.plot(independ = Color$newDepth, depend = Color$Color, type = "hist", boxp = FALSE, ylabel = "", col="gray", ylabel2 = "", xlabel = "Depth (m)")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "Brown Orange", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)
merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results##
## (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
## 0 171 42 59 61 31 10 9 2 3 1 0
## 1 50 43 53 52 30 20 9 12 16 8 5
##
## (11,12] (12,13]
## 0 1 0
## 1 0 1props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(-.22, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Color Morph Proportion")
Bars indicate relative frequency of Brown vs. Orange color morph dominance at 1m depth intervals when pooling all samples collected.
Two-Way ANOVA
Significant Interactive Effects
Dominant Symbiont by Color Morph and Reef Area
merged$Reef.Area <- ifelse(merged$Reef.Area!="Top", yes = "Slope", no = "Top")
table(merged$Dom, merged$Color.Morph, merged$Reef.Area)## , , = Slope
##
##
## Brown Orange
## C 226 122
## D 21 79
##
## , , = Top
##
##
## Brown Orange
## C 40 45
## D 12 161model=aov(Dominant~Color.Morph*Reef.Area, data = merged)
Anova(model, type = 2)## Anova Table (Type II tests)
##
## Response: Dominant
## Sum Sq Df F value Pr(>F)
## Color.Morph 21.329 1 134.208 < 2.2e-16 ***
## Reef.Area 14.489 1 91.168 < 2.2e-16 ***
## Color.Morph:Reef.Area 1.780 1 11.201 0.0008612 ***
## Residuals 111.565 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dominant Symbiont by Color Morph and Depth
Because an interactive effect of reef area and color morph was observed and slope is indicative of a depth gradient, the interaction between depth and color morph was tested here.
model1=lm(Dominant~Color.Morph*newDepth, data = merged)
Anova(model1, type = 2)## Anova Table (Type II tests)
##
## Response: Dominant
## Sum Sq Df F value Pr(>F)
## Color.Morph 24.366 1 152.751 < 2.2e-16 ***
## newDepth 9.760 1 61.187 1.909e-14 ***
## Color.Morph:newDepth 6.098 1 38.228 1.068e-09 ***
## Residuals 111.977 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Depth at which Orange switches from D to C Dominance
While brown was always C-dominated, orange was observed to switch from D to C dominance. The depth threshold at which this switch occurs was analyzed here.
threshdepth <- function(color) {
df <- subset(merged, Color.Morph==color)
plot(df$Dominant2~df$newDepth, xlab="Depth (m)", ylab = "Probability of Clade C Symbiont",
main=color)
abline(h = 0.5, lty=2)
results=glm(Dominant2~newDepth, family = "binomial", data = df)
pval <- data.frame(coef(summary(results)))$`Pr...z..`[2]
mtext(side=3, text=pval)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted ~ seq(0,11,0.01))
thresh <- ifelse(pval < 0.05,
newdata$newDepth[which(diff(sign(fitted - 0.5))!=0)], NA)
return(thresh)
}
sapply(levels(merged$Color.Morph), FUN=threshdepth)
## Brown Orange
## NA 2.75df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="sienna", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade C Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(-.22, 0), xpd = NA)
Dominant Symbiont by Depth and Reef Type
model3=aov(Dominant2~newDepth*Reef.Type, data = merged)
Anova(model3, type = 2)## Anova Table (Type II tests)
##
## Response: Dominant2
## Sum Sq Df F value Pr(>F)
## newDepth 24.789 1 123.6146 < 2.2e-16 ***
## Reef.Type 0.016 1 0.0801 0.777311
## newDepth:Reef.Type 1.650 1 8.2263 0.004253 **
## Residuals 140.774 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1df <- subset(merged, Reef.Type=="Patch")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade D Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)
Color Morph by Depth and Reef Type
model2=aov(Color~newDepth*Reef.Type, data = merged)
Anova(model2, type = 2)## Anova Table (Type II tests)
##
## Response: Color
## Sum Sq Df F value Pr(>F)
## newDepth 18.253 1 83.9552 < 2e-16 ***
## Reef.Type 0.069 1 0.3195 0.57208
## newDepth:Reef.Type 1.289 1 5.9284 0.01515 *
## Residuals 152.838 703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1df <- subset(merged, Reef.Type=="Patch")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)
Dominant Symbiont per Color Morph by Depth and Reef Type
Because depth and reef type have an interactive effect on both color morph and dominant symbiont clade, the dominant symbiont per color morph at each reef type was visualized.
df <- subset(merged, Reef.Type=="Patch")
dfo <- subset(df, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = dfo)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Patch")
dfb <- subset(df, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = dfb)
newdata <- list(newDepth=seq(0,11,0.01))
par(new = T, mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="black", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
dfo <- subset(df, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = dfo)
newdata <- list(newDepth=seq(0,11,0.01))
par(new = T, mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="red", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
dfb <- subset(df, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = dfb)
newdata <- list(newDepth=seq(0,11,0.01))
par(new = T, mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
legend("topright", legend=c("PO", "PB", "FO", "FB"), fill=c("dodgerblue3", "black", "red", "orange"), inset = c(-.2, 0), xpd = NA)
mtext(side = 2, text = "Probability of Clade D Symbiont", line = 3, cex = 1)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
Dominant Symbiont by Depth and Reef ID
model2=aov(Dominant~newDepth*Reef.Type, data = merged)
Anova(model2, type = 2)## Anova Table (Type II tests)
##
## Response: Dominant
## Sum Sq Df F value Pr(>F)
## newDepth 24.789 1 123.6146 < 2.2e-16 ***
## Reef.Type 0.016 1 0.0801 0.777311
## newDepth:Reef.Type 1.650 1 8.2263 0.004253 **
## Residuals 140.774 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1merged$Dominant2 <- ifelse(merged$Dom=="C", 0, 1)
par(mfrow=c(5,5))
domreef <- function(id) {
df <- subset(merged, Reef.ID==id)
results=glm(Dominant2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(1, 1, 3, 1))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l",
col="dodgerblue3", lwd=3, xlab="", ylab="")
mtext(side = 3, text = id)
abline(h=0.5, lty=2)
}
sapply(levels(merged$Reef.ID), FUN=domreef)
On the y-axis, a value of 1 indicates D-dominance and a value of 0 indicates C-dominance.
Mixture by Depth and Reef ID
merged$Mixture <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 1, 0)
model4=aov(Mixture~newDepth*Reef.ID, data = merged)
Anova(model4, type = 2)## Anova Table (Type II tests)
##
## Response: Mixture
## Sum Sq Df F value Pr(>F)
## newDepth 18.820 1 90.3912 < 2.2e-16 ***
## Reef.ID 9.944 24 1.9900 0.003478 **
## newDepth:Reef.ID 7.967 24 1.5944 0.036252 *
## Residuals 136.581 656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow=c(5,5))
mixreef <- function(id) {
df <- subset(merged, Reef.ID==id)
results=glm(Mixture~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(1, 1, 3, 1))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l",
col="dodgerblue3", lwd=3, xlab="", ylab="")
mtext(side = 3, text = id)
abline(h=0.5, lty=2)
}
sapply(levels(merged$Reef.ID), FUN=mixreef)
On the y-axis, a value of 1 indicates a mixture and a value of 0 indicates single symbiont clade.
Figures
Depth Influence on Dominant Symbiont and Color Morph
par(mfrow=c(3,1))
merged$DepthInt <- cut(merged$newDepth, breaks = 0:13)
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(2, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(2.1, 4, 2, 6))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="Probability of D-Dominance")
merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(3, 4, 1, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)),
xlab = "", ylab = "Probability of Orange-Dominance",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(3.1, 4, 1, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="")
df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 0, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="Depth (m)", ylab="Probabilty of D-Dominance", axisnames=FALSE, xaxs = "i", yaxs = "i")
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0, 11, 0.01), col="sienna", lwd=3)
#legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(-.13, 0), xpd = NA)
#results=table(merged$Dominant2, merged$Color.Morph)
#chisq.test(results)
#prop.table(results, margin = 2)
#par(new=T, mar=c(10, 10, .5, 6.3))
#barplot(prop.table(results, margin = 2), col = c(alpha("red", 0.35), alpha("blue", 0.35)), xlab = "", ylab = "", yaxt = 'n')
#legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.35), alpha("red", 0.35)), inset = c(0, 0), xpd = NA)Proportion of Clade D when Present
propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]
propDHist <- subset(merged, propD > 0 & propD < 1)
propDHist$propD <- cut(propDHist$propD, breaks = 5)
DCol=table(propDHist$Color.Morph, propDHist$propD)
z <- with(subset(merged, propD==0), table(Color.Morph))
o <- with(subset(merged, propD==1), table(Color.Morph))
DCol <- cbind(z, DCol, o)
par(mar=c(4, 4, 2, 0))
bars <- barplot(DCol, xlab = "Clade D Proportion", ylab = "Number of Samples", main = "", col = c(alpha("sienna", 0.55), alpha("orange", 0.55)), axisnames = FALSE, space = c(0,0.3,0,0,0,0,0.3))
axis(side = 1, at = c(1.3, 2.3, 3.3, 4.3, 5.3, 6.3), labels = c(">0.0", 0.2, 0.4, 0.6, 0.8, "<1.0"))
axis(side = 1, at = c(0.5, 7.2), labels = c("All C", "All D"), tick = FALSE)
Color Morph and Dominant Symbiont at each Reef
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["Brown.C.y"], reef["Orange.C.y"], reef["Brown.D.y"],
reef["Orange.D.y"]), radius = 7,
col = c("turquoise4", "steelblue1", "tomato", "coral"))
})
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("turquoise4", "steelblue1", "tomato", "coral"))
par(new=T, mar=c(14,21,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()
Adjusted values of Dominant Symbiont and Color Morph adjusted for the influence of depth.
Orange and Brown Color Morph
